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QJ Voronoi tessellations have been used to model the geometric arrangement 

CD of cells in morphogenetic or cancerous tissues, however so far only with flat 

hypersurfaces as cell-cell contact borders. In order to reproduce the ex- 
perimentally observed piecewise spherical boundary shapes, we develop a 
r~| consistent theoretical framework of multiplicatively weighted distance func- 

Q-i tions, defining generalized finite Voronoi neighborhoods around cell bodies 

of varying radius, which serve as heterogeneous generators of the resulting 
model tissue. The interactions between cells are represented by adhesive and 
^ repelling force densities on the cell contact borders. In addition, protrusive 

locomotion forces are implemented along the cell boundaries at the tissue 
margin, and stochastic perturbations allow for non-deterministic motility ef- 
fects. Simulations of the emerging system of stochastic differential equations 
Q_i for position and velocity of cell centers show the feasibility of this Voronoi 

method generating realistic cell shapes. In the limiting case of a single cell 
CN pair in brief contact, the dynamical nonlinear Ornstein-Uhlenbeck process is 

analytically investigated. In general, topologically distinct tissue conforma- 
\^ tions are observed, exhibiting stability on different time scales, and tissue 

coherence is quantified by suitable characteristics. Finally, an argument is 
derived pointing to a tradeoff in natural tissues between cell size heterogene- 
ity and the extension of cellular lamellae. 
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1 Introduction 

A Voronoi tessellation is a partition of space according to certain neighborhood relations 
of a given set of generators (points) in this space. Initially proposed by Dirichlet for 
special cases |16| . the method was established by Voronoi more than 100 years ago |45j . 
The geometric dual of the Voronoi tessellation was proposed by Delaunay in 1934 - 
and therefore is called Delaunay triangulation. It connects those points of a Voronoi 
tessellation that share a common border. Since the latter can be directly constructed 
out of the former, both terms are sometimes used equivalently. In the following years, the 
method was rediscovered throughout other fields, which accounts for many other names 
designating the very concept, such as Thiessen polygons [43] in meteorology or Wigner- 
Seitz cells [48] and Brillouin zones [13] in solid state physics. With the technological 
and scientific advance, the method became feasible in computational geometry [39] , and 
since then has widely evolved, cf. (9J, [8] , making it appealing for biological applications. 

In particular, Voronoi tessellations have been applied to represent various aggregates of 
cells and swarming animals. Initially, Honda proposed the method in two spatial dimen- 
sions [24J. The first applications to biological tissue were cell sorting simulations, however 
starting from artificially shaped quadratic cells [41]. Then, morphogenesis and its un- 
derlying intercellular mechanisms were studied starting from a pure Delaunay mesh and 
simulating vertex dynamics |47} 146] , yet without using Voronoi tessellations explicitly. 
In contrast, by applying transformation rules like mitosis combined with Monte-Carlo 
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dynamics, evolving multicellular tissue was represented by Voronoi tessellations [18]. In 
particular, growth instabilities, blastula formation and gastrulation could be conceived 
within this framework [T7]. Similar effects were reproduced by using vertex dynamics 
|14j . Moreover, cell organization in the intestinal crypt was modelled using spring forces 
and restricting the motion to a cylindrical surface [32 . An application to bird swarming 
together with the proposal of a continuum formulation was given in [3|. Finally, the 
influence of shear stress on the evolution of two-dimensional tissues was studied p~5] . 

Only quite recently Voronoi tessellations have been extended to be used as a model for 
three dimensional tissue, again using vertex dynamics [25 . Other authors use optimized 
kinetic algorithms |36| [TTj to employ generalized Voronoi tessellations (discussed as dif- 
ference method in this article), with cell-cell and cell-matrix adhesion |37j . Marginal cells 
have been closed by prescribing a maximal cell radius, enabling the study of the growth 
dynamics of epithelial cell populations |21| . So far, however, the cell-cell boundaries 
were exclusively represented by flat hypersurfaces. 

In contrast, when observing two-dimensional monolayers of keratinocytes, for example 
cf. |49| \31\ |4"4"] . the cell-cell contact borders visible from staining cadherin-complexes 
frequently appear as circular arcs, whose shape and length seems to be determined by 
the constellation and size distribution of neighboring cells. Moreover, the forces between 
such cells are influenced by filament networks or bundles meeting at the cell-cell junctions 
and eventually balanced by elastic counterforces [U IH [291 SO]- Therefore, a geometrical 
and dynamical modeling framework is required that reproduces the observed cell shapes 
and simultaneously allows for quantifying the cell-cell interaction forces as well as the 
active locomotion forces appearing at the free cell boundaries. Here we present a simple 
and effective solution of this task by using a suitably weighted Voronoi tessellation. 

This article is organized as follows: In section [2] Voronoi tessellations are introduced 
in a general manner. Next, two types of weighted square distance functions are in- 
troduced, using the method of difference and quotient, respectively, and their particular 
consequences for cell tissue modeling are investigated. Inspired by the intricate interplay 
between cytoskeletal filament bundles and cadherin-catenin cohesion or integrin adhe- 
sion sites, the forces on the intercellular and exterior cell borders are proposed in section 
[3] after discussing the emergence of cell shape within our model. Then the dynamics of a 
whole cell aggregate is defined, directly leading to analytical results on cell pair contacts 
in section |4j After simulation studies of meta-stable states during tissue equilibration 
and robustness of tissue formation under the influence of various model parameters in 
section [5] we conclude with a discussion of our results in section [6j 

2 Generalized Voronoi tessellations 

Let {gi : i = 1 . . . N} denote a finite set of N generators or points Xj in Euclidean, 
n-dimensional space W 1 . 

Definition 1. The Voronoi cell of a generator gi = x, is defined as 

Vi = {x G R n : Pi(x) < P,(x) Vj/i}, (1) 
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where V% (i = 1 . . . N) is a given set of continuous, generalized square distance functions 
on W 1 with the property that \/i : Xj £ Vj. 

Thus, Vj represents an open neighborhood of Xj, containing all points x that are 
"P-closer to Xj than to any other xj. 

Definition 2. The contact border between two points Xj and Xj is defined as the inter- 
section of the closures of Vi , Vj : 

Tij = Vi n Vj with i j. (2) 

The total boundary of the Voronoi neighborhood around Xj then is 

dVi = \jT tJ . 

The contact border therefore is the set of all points P-equidistant from Xj and x.,-, 
namely 

r ij ={xef 1 :P,(x)=P ] (x)<? t (x)V^y}. (3) 

The Voronoi tessellation in general form is then given by {Vi,Tij; i,j = 1, . . . ,N} and 
covers the whole space, where so-called marginal neighborhoods Vi extend to infinity. 
Depending on the particular choice of the generalized square distance functions, Tij can 
take various shapes. In the standard Euclidean case, 'Pj(x) = |x — Xj| 2 , the contact 
border Ty is the perpendicular bisector of the line segment from Xj to x^, an (n — 1)- 
plane. Then Vi is bounded by a convex, not necessarily finite polytope and is called the 
(classical) Voronoi neighborhood |45j or Dirichlet domain [16 . 

The modeling aim here is to represent biological eukaryotic cells in connected 3- 
dimensional tissues or confluent 2-dimensional cell monolayers as Voronoi neighborhoods 
Vj, as in the 2-dimensional pictures in figure [TJ In a minimal approach we define the 
points Xj as centers of the visible, mostly ball-shaped cell bodies, containing the cell 
nuclei plus other cell organelles such as mitochondria or the Golgi apparatus. By at- 
tributing a finite radius rj > to each Xj, the Voronoi concept is extended to generators 
gi = B n (xi) of positive finite volume, being a suitable representation of cell bodies. Since 
these are rather solid in comparison to the rest of the cell, it is assumed that the B ri (xi) 
do not overlap. Then Vi \ B ri (xi) represents the protoplasmic region of the cell i, which 
for n = 2 appears as a fiat lamella in light microscopy. Clearly, the natural condition 
£>r;(xi) C Vi requires that 

Vx€^(x7) Vj^i: ^(x)>^(x). (4) 

It shall be seen later, that this condition is fulfilled for the chosen generalized square 
distance functions. 

Furthermore, certain weights Wi G M + , on which the square distance function Vi may 
depend, are assigned to each cell i. These weights Wi = w{ri) are assumed to be strictly 
monotonically increasing functions of the cell body radius r^, with the intended effect 
that stronger weights induce larger cell sizes, by shifting the Voronoi contact borders 



4 





Figure 1: (a) Micrograph of human keratinocytes in a section of stratum spinosum in 
vivo (reproduced from [19]), and (c) phase contrast microscopic photograph of 
human keratinocytes in-vitro (by courtesy of Institute of Cell Biology, Bonn 
University). Image (c) is extracted from the supplementary movie movO . avij of 
the dynamics at the edge of an almost confluent monolayer in a wound scratch 
assay. White arrows in (c) indicate a round cell (upper) and a cell pair com- 
peting for its influence region (lower). According to these observations, tissue 
cells in a two-dimensional geometry are modelled as Voronoi neighborhoods Vi 
containing the ball shaped cell bodies 23 r .(xj) surrounded by a so-called lamella 
(b). 



outwards. Importantly, different choices of how Vi depends on it;, could lead to different 
cell shapes. Out of many possible generalizations of Voronoi tessellations (for a review 
see [H]), we only discuss two straight-forward ways here, which are determined by the 
set of all cell center positions, body radii and weights {gt = (B ri (xi), Wi)}. 

2.1 Difference method 

The partition of space into cells is obtained by a Voronoi tessellation using the Euclidean 
square distance function with subtracted weights 

7>i(x) = (x- Xi ) 2 -™ 2 , ( 5) 

which has previously been used in |24[ |23] without and in |37| [TU] with weights Wi. 

For the following we denote x™ = (x, + x,-)/2 the cell pair mid-point, dij = |xj — Xj\ 
the Euclidean cell center distance, and d,j = (x, — Xj)/dij the unit vector of the oriented 
axis connecting the two cell centers. Thus, from equation ([3]) the condition for a point 
x to be located on the contact border Tij reads as 

— 

(x-x ij )-d i j = — ( 6 ) 

being equivalent to a linear hyper-plane equation. As for the classical Voronoi partition, 
the contact (n — l)-plane between two neighboring cells i, j is perpendicular to the vector 
connecting the cell centers. However, now the position of the contact border plane along 
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the connecting vector, and thereby the sizes of the Voronoi cells, depend on the weights. 
In figure [2] the geometry of a cell pair with its separating border is illustrated. 




Figure 2: Geometry of a cell pair and its border using the difference method. There is a 
distinct direction given by the axis line connecting x,- and Xj. The circular 
cell bodies around Xj, Xj and their radii r^, tj are indicated by green lines, the 
contact hyperplane by a red one. Other quantities are explained in the 
text. 

In order for a distance function to yield a partition describing an aggregate of biological 
cells in living tissue, the border Tij has to be located between the surfaces of the non- 
overlapping cell bodies. With the corresponding distances as denoted in figure |2j this is 
equivalent to 

Si > A 6j> 0. (7) 
These constraints have consequences for possible choices of the weights: 

Lemma 1. Let {Vi,Tij} be a Voronoi tessellation of non- overlapping generators {gi = 
(B ri (xj), Wi)} constructed from V% according to the difference method in definition |5p 
with positive weights Wi. Then the inequalities |?p are satisfied for all cell pairs i j 
with arbitrarily small but positive cell body distance 5ij = 5i + 5j, if and only if 

Vi : Wi = n. (8) 

Proof. : The contact border equation V\ = Vj evaluated at the intersection point x = 
Sij, see figure [2J yields (rj + 5i) 2 — wf = (rj + 5j) 2 — Wj. Together with Si + Sj = 5ij we 
obtain the representation 

<5 ij (^+2r J ) + K 2 - U ; 2 )-(r2-r 2 ) 
- — > 0. (9) 



2(n + rj + 8^) 

Since 5ij > can be arbitrarily small for fixed ri,rj and Wi,Wj, the condition <5j > 
implies wf — rf > — r 2 . By exchanging i and j, the second condition 5j > enforces 
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the equality and the existence of a joint constant C with 

W k = r l + C for k = h3- 

Since the definition [T] of a Voronoi cell is independent of such an additive constant in 
equation §5§ , we can set C = and obtain the result of the lemma. □ 

Thus, while a Voronoi tessellation can be formally defined using arbitrary subtrac- 
tive weights, the constraint of non-overlapping cell bodies leads to the unique choice of 
weights Wi = Ti. Furthermore, these weights imply a simple characterization of the cell 
bodies £> ri (xj) = {x : 'Pj(x) < 0}, so that inequality Q is fulfilled under the assumption 
B n (xj) C\B r . (xj) = 0. An illustration of a two-dimensional Voronoi tessellation with such 
weights is shown in figure [3^a). The geometric interpretation of this choice of weights 
gives rise to the "empty orthosphere criterion" for a regular triangulation in \37\ [35] [TO] , 
since the squared radius of the "orthosphere" equals the P-distance of three or more 
neighboring generators from their common Voronoi border junction, consisting of red 
lines in figure [3^a), compare the analogous figure 1 in [TO] . 

From equation ([9| and condition Q we obtain the dependence 

= (%(jg_+jrj) 5i_ = 5jj + 2rj 

2{ri + rj +Sij) 5j 5ij+2ri' 

For one, if > rj (as in figure [2]), then <5, < 8j. Furthermore, for fixed 5{j > and rj, Si 
is monotonically decreasing in r{. This means that for growing cell body radius r» > rj, 
the distance Si between cell body and contact border Tij (attained at Sy) would shrink, 
thus also the size of the protoplasmic region Vj \ S ri (xj). However, such a behavior 
is contradictory to empirical observations: If two cells i and j touch each other, then 
the cell with a larger cell body should also have a wider cytoplasmic region along the 
contact border, see [2] (figure 5E) and figure [TJ Therefore, the difference method is not 
appropriate and an alternative method is required. 

2.2 Quotient method 

In the previous section it was found that, with subtractive weights in the "P-distance of 
the generalized Voronoi tessellation, the emerging cell contact border are planar surfaces. 
In contrast, if one divides the Euclidean distance by weights, then the cell contacts are 
spherical with the generalized square distance function defined as 

^(x) = (10) 

Having its roots in computational geometry (see [7] and references therein), this method 
was introduced as a model for attraction domains of restaurants |6i more than 20 years 
ago. To our knowledge, it so far has not been used for physical or biological applications. 

For simplicity of calculation, let the midpoint Xjj := be the origin of the coordinate 
system, while dfj remains the oriented cell-cell axis (see figure kl). Starting from 'Pi(x) = 
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Figure 3: Comparison of two generalized Voronoi tessellations from the same set of gen- 
erators. While the difference method (a) yields polygonal cells the quotient 
method (b) yields cells with piecewise spherical boundaries. In both cases, the 
weights are given by Wi = r{. The cells are described by their cell center (green 
star) and their body (thick green circles). The Voronoi borders between cells 
are red lines, while the neighbor relations (i.e. those cells that share a common 
border) are indicated by thin green lines connecting their centers. 



"Pj(x) with Wi ^ Wj one arrives at the equivalent condition 

Uj) = Rij 



(x-M«) 2 = i& (11) 



for the point x to be on the border Tij. Clearly, equation (11) describes an n-sphere 
around 

wf + W 2 a WiWj 

M« = q iyXi with radius Ra = - — * ozda, (12) 

K-w) ("7 w j) 

where dij = |xj — Xj| = 2|xj|, resulting in a so called circular Voronoi tessellation. 
Assuming n > Tj (as in figure |4]), also the weights fulfill Wi > Wj according to our 
monotonicity assumption on Wi = iw(rj). Thus, from equation (12) the center My of 
the hypersphere Ty is always situated on the side of the cell j with the smaller radius 
Tj. The contact sphere intersects the cell center connection segment Xj , Xj at a unique 
contact point determined by 

Sy-dy "•' "•'•' / ". (13) 
3 3 Wi + Wj 2 v ' 

Similar as for the difference method (see figure |2j, Sy is situated on the side of the 
smaller cell body from the mid point Xy = 0. Thus, the contact sphere contains the 
body of the cell with smaller weight, as indicated in figure |4j Once both weights are 
equal, equation ^ and thereby equation (13) simplify to 

x • dy = Sy • dij = for Wi = Wj, (14) 

revealing Ty as the classical Voronoi bisector line without weights. 

In analogy to Lemma [T] for the difference method, we obtain the same unique specifi- 
cation of weight functions here as well: 
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Figure 4: Geometry of cell pair and its pair contact border using the quotient method. 

In contrast to the difference method, this contact border Yij is a sphere around 
Mj,- with radius R;j. Its two-dimensional section is drawn in a red line. 



Lemma 2. Let {Vi,Tij} be a Voronoi tessellation of non-overlapping generators {gt 



(B ri (xi),Wi)} constructed from V% according to the quotient method in definition (10) 
with positive weights W{. Then the inequalities |?j) are satisfied for all cell pairs i ^ j 
with arbitrarily small but positive cell body distance Sij = Si + Sj if and only if 

Vi : Wi = n. (15) 

Proof. : The contact border equation V% = Vj evaluated for the point 

figure |4j yields (rj + 5i)/wi = (rj + 8j)/wj. Together with 5i + Sj = Sij we obtain the 

representation 

Wi rjwi - nwj 

Si = Sij ; ^~ ; ~>0. (16) 

Wi + Wj Wi + Wj 

Since Sij > can be arbitrarily small for fixed r(,Tj and W{,Wj, the condition Si > 
implies riWj > rjWi. By exchanging i and j, the second condition Sj > enforces 
equality and the existence of a joint positive constant with 

Wk = r k -C for k = i,j. 

Since the definition [T] of a Voronoi cell is independent of such a multiplicative constant 



in equation (10), we can set C = 1 and obtain the result of the lemma. □ 



Thus, further on we can choose the weights Wi = rj when using the quotient method. 
Then the cell bodies are characterized as B n (xj) = {x : Vi{x) < 1}, so that inequality 



4| is fulfilled under the assumption B ri (xi)PiB rj (xj) = 0. The emerging circular Voronoi 
tessellation is illustrated in figure [3^b) . Rewriting equation (16) for i and j yields 



n + rj 
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and thus the regular partition property 



n 



(17) 



meaning that the partition of the distance between cell bodies £> r .(xj) and B r -(xj) by 
the contact arc is proportional to the ratio of body radii. Most importantly, in contrast 



to the difference method, the "P-distance in equation (10) ensures that Si monotonically 
increases with r^, while 5j = 5a — <5j decreases, if r,- > and Sij > are held fixed. As 



a consequence, the cell i grows with increasing r^. This property can also be observed 
for in- vivo cell monolayers, see figures [TJ [5] and figure 5E in |44| . Thus, for the further 




Figure 5: Typical microscopic picture of an epithelial monolayer (left) cultured of human 
keratinocytes, cell nuclei (green, black) and cell-cell contacts (yellow, red) are 
visualized by suitable staining (reproduced from [31] ). In the simulated cell 
tissue (right) the sizes and positions of cell bodies (green) have been roughly 
adapted to deliver the observed contact arcs (red). 

discussions in this paper, the quotient method will be used exclusively to define Voronoi 
cells. Yet the tessellation is still unbounded due to infiniteness of marginal cells. 



2.3 Closure of the Voronoi tessellation 

In order to avoid infinitely extended cells at the tissue margin, the initial definition [T] 
of a Voronoi cell has to be modified. To this end, the method of finite closure for the 
difference method, as used by Drasdo and coworkers [181 EI] ; 1S extended to be generally 
applicable. 

Definition 3. Let P max £ M + . The finite Voronoi cell of a generator (B ri (x.i),Wi > 0) 
is defined as 

V, = {xGM n : 7>i(x) < min (7^(x), P max ) Vj + i) . (18) 
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The exterior boundary closing a marginal Vi is 

T i0 = {x G M n : n(x) = 7W \ U Vj, (19) 

and the total boundary of the Voronoi neighborhood around x, is 

dVi=T i0 u[jT ij , 

where now the contact border between cell i and j is given by 

Tij = {x6f: Vi(pt) = Pj(x) < min(n-(x),P max ) VA; + i,j}. (20) 

For any choice of "P max > 0, the Voronoi tessellation generated from a finite set of cell 
bodies {(B ri (x.i), u>i)} comprises a bounded region of the whole space, representing a cell 
tissue of finite extension. In figure [6] we depict the model representation of an in vivo 




-5 5 10 15 20 



Figure 6: Model representation of a bounded cell monolayer using the quotient method 
with "Pmax-cutoff. Axis tics are in units of /xm. 

cell monolayer (see figure [5| generated by the quotient method using definition (10) with 
Wi = T{. The exterior boundaries r,o of the Voronoi neighborhoods for marginal cells 
are circular arcs drawn as black lines. Apparently the size of a cell i is influenced by 
the two parameters r% and P max - While P max regulates the overall cell size by globally 
scaling each rj, the ratios ri/rj determine the partition of space between each cell pair 
i,j by specifying the actual position of r^-. We remark that both ? max and {rt} are 
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Figure 7: Maximal P-distance P max and its effect on the neighbor relation of two cells 
The cells i,j can only share a common T{j and thus be neighbors, if there 
is a non-empty overlap region (shaded) given by the intersection of their free 
balls Bf^Xj), ,Br. (xj) (blue outer circles). The other quantities are explained 
in the text. 



accessible to experimental determination using image analysis tools, see figure [5] and 
especially figure [TJc). 

The necessary condition Pi(x) < P max for a point x to be within cell i defines a ball 
BR i0 (x-i) around Xj, which will be called free ball further on. It has the squared radius 
i?f = wf + P max in the difference method and Rf = wf • P max in the quotient method. 
As illustrated in figure [7j this has the effect that two cells i,j can only be neighbors, if 
their free balls S/j i0 (xj),B^ (xj) overlap. In case of no contact, Vi = i3ft i0 (xj) represents 
an isolated spherical cell, which clearly has to include its cell body i3 n (xj). Thus, the 
condition 

9 o I r f + ^max (difference method) / v 

r? < R% = I 1 (21) 
I r i • "Pmax (quotient method) 



is imposed, meaning that Pmax has to be chosen so that 

Pmax > (difference method) 
Pmax > 1 (quotient method). 



(22) 



In this way, the larger P maX ) the larger is the radius of isolated cells relative to their cell 
body. 

Within the difference method, this P max -closure is straight-forward because the planar 
cell contacts lead to starlike (even convex) Voronoi cells Vi. Recall the definition of 
starlikeness with respect to the center x«: Vx G Vi also x^x C V,. In fact, a P m ax- 
closed generalized Voronoi tessellation has been applied to epithelial tissue modeling by 
prescribing Z3ij i0 (xj) for each cell [21]. However, within the quotient method the situation 
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is more complicated. In analogy to the difference method it is reasonable to require that 
the Voronoi cells Vi are star-like domains with respect to Xj. In order to ensure star like 
cells within the quotient method, ? max must not be chosen too large. Consider the cell 
pair as sketched in figure [7j The straight line connecting Xj and Ct is a tangent to Tjj. 
Thus it is clear from the geometry that both Vi and Vj are star-like domains with respect 
to Xj,Xj, if their corresponding free balls S^ i0 (xj), Z3p (x,) do not extend beyond the 
point Ct- Before we proceed, we introduce the cell size homogeneity quotient 

iy* . _l_ rf* . rf* _l_ iy* 

s*, ■ 1 1 T ' 7 ' max ~r ' mm /no\ 

Q = mm — = , (23) 

h3 | f% Tj | ''max ?"min 

where the last equality follows from monotonicity arguments. Therefore, Q = Q({ri : 
i = 1 . . . N}) is a measure of the uniformity of cell sizes within a tissue, with Q = oo for 
equal r» and Q as 1 for r max >> r m i n - 

Proposition 1 (Starlike cells). For a finite Voronoi tessellation generated from non- 



overlapping {£> r .(xj)} by using the quotient method in definition (10) with weights Wi = 
r{, the resulting Voronoi cells Vi are starlike with respect to Xj, if the maximal V- distance 
"Pmax fulfills the homogeneity constraint 

1 < Pmax < Q. (24) 

This condition on the tissue properties will be crucial later on and guarantees that 
each actin fiber bundle emanating radially from dB ri (x.i) intersects the boundary dVi 
only once, see figure M 



Proof. : From fundamental trigonometric relations follows the angle Z(T, x,-, Ct), namely 

(25) 



b T = - 
1 2' 



and geometric similarity of the triangles A(Mjj,Xj, Ct), A(xj, Ct, xj). Thus, with Vj < 
Ti, we have for point A = Ct 

cos 9 T = — (26) 

n 

RjT = j - ■ dij, RiT = , - ■ dij. (27) 

rf-r]\ \J\ r i- r j\ 

With the last two equations, the maximal distances of a point A on from the cell 
centers have been identified for each cell pair. Starlikeness of Vi is equivalent to the 
condition R 2 i0 < R] T , where R? = V max rf, so that P m ax < df-/\rf - r||, which can be 



fulfilled by requiring T^max — Qi since \/i,j : (n + rj) 2 < dfj. With the condition fc2j 
the assertion follows. IJ 

In particular, starlikeness prohibits engulfment of one cell by the other, so that Bn i0 (xi) 
may not contain Bn j0 (xj) completely for r{ > rj. Note that within sufficiently large tis- 
sues, the smallest and biggest cell will usually not be in contact, which relaxes inequality 
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Figure 8: Electron microscopic picture of a keratinocyte (courtesy of Gregor Wenzel). 

There are two distinct regions within the cell. An inner, more dense region 
with ruffled membrane structures appearing in white (cell body), and an outer, 
flat region with larger ruffles or filopodia near the cell margin (lamella), which 
in the case of no contact with other cells forms a ring of more or less constant 
diameter around the cell body. See also figure [T] 



(24) into the condition: 



KV max < min ri + r i = : Q nh . (28) 

neighbors i,j | Tj — Vj \ 

In such a circular Voronoi diagram there may be 0(N 2 ) cell-cell contacts and vertices, 
in particular for low cell size homogeneity quotients Q. The supplementary material 
contains the extensively commented GNU octave routine mwvoro.m (Matlab compati- 
ble) used to create and visualize a circular closed Voronoi tessellation, together with 
configuration and plotting facilities. The partition of space into distinct cells and their 
neighborhood relations has an algorithmic complexity of 0(N 2 ), N being the number of 
cells. Due to the extension of closing marginal cells by 'Pmax-arcs, we do not follow the 
method proposed by Aurenhammer and Edelsbrunner [7 ]. In particular, we do not use 
polyedral "cell complexes" in an inverted three-dimensional embedding of the Voronoi 
generators. Nevertheless we retain the same optimal algorithmic complexity. 



3 Cell shape and dynamics 

When a single cell is placed on a two-dimensional and adhesive substratum, it usually 
spreads into all directions attaining a circular shape like a fried egg, as can be observed in 
figure [8] Thereby, the exterior visco-elastic lamella along the free boundary Tjo, which 
consists of parts of £>R i0 (xj), supports the protruding and retracting cell edge. This 
smooth flat region contains a network of dynamic actin filaments situated around the 



14 



inner, almost solid cell body B ri (xj), see e.g. [U 20]. The maximal spreading radius 
Rio = V^max • r « is determined by the strength of adhesion to the substrate and the 
equilibrium between protrusive activity at the cell periphery I^o and the contractile 
retrograde actin flow [30 . Assuming that for given adhesiveness the averaged local 
cytoskeletal network volume fraction (9 in [30]) in the lamellae has a certain value go 
independent of cell size, the weighting of proportional to the cell body radius follows. 

Usually, the irregular activity of living cells at their lamella edge leads to a curled 
cell boundary, see figures [ijc) and [8] Neglecting fluctuations on short time O(10s) and 
length scales 0(l/mi), the free portions of the cell edge I^o are approximately taken as 
circular in this model. The actual cell edge fluctuates within the vicinity of the smooth 
arcs, representing the averaged position of the plasma membrane, and will later on be 



taken as the source of stochastic perturbation forces, see section 3.4 Thus, in our model 



the active lamella region of a single free cell is approximately ring shaped and has a 
width of 

S i0 = (Vr^c - l)r<. (29) 

Once two epithelial cells i, j come close enough to interact, the two adjoining lamel- 
lae compete for the occupation of the region in between them. Eventually they form a 
contact border, which exhibits microscopic fluctuations due to local plasma membrane 
flickering. Yet it approximately attains the shape of a circular arc, whereby small gaps 
between the two cell membranes are neglected. This experimental fact (see e.g. corre- 
sponding figures in |49| |3"T] |H] and [5] [T] in this article) is well represented by the Voronoi 
border Tij resulting from the quotient method defined by equation (10). In this way, 
within our tissue model, the cell boundaries are merely composed of piecewise circular 
arcs, and the cell bodies are not necessarily located in the middle of the cells. By suit- 
able choice of Wi (Lemma [2| there is always some lamella region separating the cell body 
from the neighbor cell (Si > 0, also cf. figure [4]). 



3.1 Interaction forces between cells 

The cytoskeleton with its network of filaments often features bundled structures, which 
are commonly visible as so-called stress fibers, emanating from the cell body or nucleus in 
radial direction. According to [I], bundles of filamentous actin attach to transmembrane 
complexes called adherens junctions, which are made from e.g. catenins on the cytosolic 
side and cadherins at the exterior of the cell. Furthermore, intermediate filaments such 
as the rope-like keratin tie in with rivet-like desmosomes at the cell membrane. By 
connecting neighboring cells, these structures stiffen and strengthen the tissue coherence. 
For example, in certain epithelia, cadherin-catenin adherens junctions comprise a whole 
transcellular adhesion belt. 

Inspired by this observation, it is assumed that the attractive force between two cell 
bodies B ri (xj),B r .(xj) is transduced by radial filament structures extending towards 
the cell boundaries. Thereby, the filaments of one cell connect to those of the other 
and form pairs along the contact border Tij. Thus, the intercellular adherens junctions 
emerge according to the respective filament densities as emanating from cell i and j, 
respectively (see figure pj. Furthermore, the connecting cell-cell junctions are not fixed 
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Figure 9: Pairing of filaments from one cell to the other cell. For the definition of symbols 
and angles to describe the geometry of pairing filaments of a neighboring cell 
pair see figure [7j 

but undergo dissociation, diffusion, and renewed association. Motivated by protein (e.g. 
cadherin) diffusion properties in membranes \22\ [12], this process is considered to be 
fast (seconds) compared to the slower time scale (several minutes) of cell deformation 
and translocation. In this way, pair formation of cross-attachments between filament 
bundles from both cell bodies can be regarded as a pseudo-stationary stochastic process 
|19j . In order to compute the interaction force between two cells, one needs a suitable 
expression for the density of pairing filaments p{6) at the border of the cells 

3.2 Filament pair density at cell-cell contacts 

Consider the cell pair as illustrated in figures [7| an d [9] with cell body radii r{ > rj and the 
distinguished Voronoi weights as in equation ( |15| . Let 6 parameterize the contact arc 
given by My , Rij , and let A be the corresponding point upon that arc. Starting from 
the surface of the cell bodies B ri (x.i) and B rj (xj), filaments extend in radial direction 
under angles fa(0) and <j>j(0), respectively, to eventually meet at A. Furthermore, let Ri 
and Rj denote the distances between the cell centers and A. The density of filaments is 
assumed to be constant on the surface of each cell body, given by a universal value p > 0. 
In order to construct the pairing density of filaments p{0) along the contact surface, these 
cell body surface densities are mapped onto by equating the corresponding surface 
elements 

Pi(d)Rijd9 = prid(f)i, pj{9)Rijd6 = prjdcpj. (30) 
With A = A(0) G r^, (cf. figures and it holds 

Rij sin 6 = Ri sin fa = Rj sin cpj (31) 
Rij cos = \xj — Mjj | + Rj cos <pj. (32) 
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The defining condition for the contact border in equation ^ can be written as 



Rj = rjRi 



with T) 



W: 



W; 



< 1. 



r, 



Thus, from equation (31) we obtain the simple relation 

sin (j>i = 77 • sin <j)j 



(33) 



(34) 



between the two angles (f>i(9) and (f>j(9), so that differentiation with respect to 9 yields 
the proportionality 



COS ( 



i] 



COS ( 



d9 



d9 ' 



(35) 



where n r , 



.) = Vl + (l-r/ 2 ) 



tarn 



Moreover, by solving equation (31) for Rj in 



terms of Rij, inserting it into equation (32), and using the relations (12) we get an 
explicit expression for tan <j>j in terms of 9 



tan ( 



sm a 

cos — rj ' 



(36) 



which holds for all \6\ < Ot, with cos 9t = Tj, or equivalently, \<pj\ < vr/2, see equations 
( 25p6 ) . Finally, by differentiation of equation ( 36 ) with respect to 9 we obtain 



tan ( 



d6> 



1 + tan 2 , 



tan + 



1 



1 — 77 cos 9 
tan 9 J 1 — 2r/ cos 9 + r\ l 



> 0. 



(37) 



It is assumed, that the pairing density function p(9), depending on pi{9) and Pj(9), 
is even in 9, maximal at 9 = 0, and strictly monotonically decreasing for increasing \9\. 
Here, two exemplary models to specify such a density function p(ff) are discussed: 



Model 1: Minimal density pairing. If locally one cell has less filaments binding to T^- 
than the other, then there will be a pairing match for all of its filaments. Thus, the local 
density of pairs on Tjj will equal the lower filament density: 



mm 



R 



p . / dcpi d(j)j 



'j 



(38) 



where we used the identities (30). With k, > we conclude from equation (35) that 



p{9) = Pi{9) < Pj{9). Therefore, an explicit representation of p in terms of <j>j and its 
derivative is 



Pmm(9) = Pi (9) = p- 



R 



1 1 " K rj(<fijj 



d9 



The emerging behavior of p(9) is visualized in figure [LOj where the maximal angle 9t 



as found before in equation (26), can be clearly seen in plot (b). It does not appear so 



prominent in plot (a), where it varies with n 
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Figure 10: Plot of density function p(9), described by model 1 (minimal density pairing). 

(a): The distance of cell bodies 5{j > and Tj > are fixed constants, while 
r-i > Tj successively increases, (b): rj > r j > are fixed constants, while the 
distance of cell bodies successively increases. 



Model 2: Mean density pairing. Assuming that each filament from either of the neigh- 
boring cells has a probability to randomly form a pair at some junction on Ty, the 
resulting pairing density can be defined as the geometric mean of pi and pf 



(e) = Jp i (6)- Pj (e) = p 



Rij ■ V^-vi^j) d0 



(39) 



In figure 11 the emerging behavior of p(0) is visualized, showing an even more expressed 
cut-off at 9 = 0t 
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Figure 11: Plot of density function p(0), described by model 2 (mean density pairing). 

(a): The distance of cell bodies <5jj > and r,- > are fixed constants, while 
ri successively increases, (b): ri > and rj > are fixed constants, while the 
distance of cell bodies successively increases. 



From figures 10 and 11 it becomes apparent that p(9) is even in 9, maximal for 9 = 0, 
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and strictly decreasing for increasing \9\ in both methods. Whereas model 1 captures 
the maximal filament pairing that could be realized for long term association at fixed 
adherens junctions, model 2 describes the pseudo-steady state of short term stochastic 
filament association, which will be considered further on. 



3.3 Pair interaction force 

Consider two cells i, j touching each other, so that if their cell body distance 5™ = 5i + 5 



was further increased, they would dissociate. By equation (29) this limiting rupture 

ij 



distance <5,-™ p ^ is given as 



^ UP) = Sio + S j0 = ( - 1) (n + r 3 ). (40) 

Then according to the previous assumptions, any paired couple of actin fibers meeting 
at an adherens junction in the contact boundary Yij near the intersection point 
(see figure |4| develops a certain positive stress between the two cell bodies. According 
to the assumption made at the beginning of this chapter, this stress depends on the 
mean volume fraction qo of the contractile cytoskeletal network, which before touching 
was equal in both contacting lamellae of width Sio,Sjo, respectively. If now Sy further 
decreases, then both lamellae will be compressed by the equal factor Si/Sio = Sj/5jQ = 



5ij/5ij Up ^ < 1 as a consequence of the Voronoi partition laws (18) and (19). Thus, we 



can suppose that (a) the mean volume fraction in both lamellae increases to the same 
value q satisfying the inverse relation 



^(rup) 



q 5. 



(41) 



and (b) any paired actin fibers develop the same stress between their adherens junction 
and the corresponding cell body, with a strength / = f{q) that, for simplicity, depends 
only on the common cytoskeletal volume fraction q. Since the cytoskeletal network 
consists not only of cross-linked actin-myosin filaments but also of more or less flexible 
microtubuli and intermediate filaments (as keratin, for example) |29 | 142 1 I4"0" l IT], the stress 
function f{q) has to decrease to (large) negative values for increasing q — > (/max — 1 • Here 
we chose the simple, thermodynamically compatible strictly decreasing model function 

/(<?) = /int (ln(l - q) - lng - In z c ) . 
The corresponding convex generalized free energy T satisfies 

^1-9) = (1 -?)(/(?) "/int) 

for < q < 1 (cf. [2]), where the positive constant z c < 1/qo — 1 determines the critical 



volume fraction q c = 1/(1 + z c ) > qo such that f(q c ) = 0. Applying transformation (41 ) 
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we finally obtain an actin fiber stress function that depends only on the relative cell 
body distance Ay = <%/£y Up ^ < 1, namely 



/(A^) = / int • In 



Ay - A r 



A 



crit 



A r 



(42) 



where < A min = q < q {l + z c ) = A crit < 1. 

The derivation of this stress model relies on the simplifying assumption that according 
to equation (41 ) the stress of each paired filament extending from cell body B ri (xi) to the 
adherens junction at Ty is completely determined by the adhesion strength (appearing 
as coefficient /i n t) and the cytoskeletal state q of the intermediate lamella near the 



horizontal cell-cell connection axis in direction dy, see figures |7jand|9| Moreover, relative 
to this coordinate frame the paired filament orientations are Rj = (— cos</>j,sin<^j) and 
Rj = (cos 4>j, sin^)j), so that the corresponding adherens junction at Ty experiences 



two force vectors f« = —/(Ay) • Rj and fj 



-f(Aij) ■ Rj with opposing horizontal 



components. However, their resultant vector fj + fj generally does not vanish (except 
for fa = <j)j = 0). It has a negative vertical component —/(Ay) • (sin^j + sin cpj), which 
could pull the adherens junction towards the cell-cell connection line along the contact 
boundary Ty. 

Therefore, some counterforces due to substrate adhesion via e.g. integrin |20| [23] or 
frictional drag have to be supposed in order to guarantee the assumed pseudo-stationary 
equilibrium condition for Ty. Using the simplifying decomposition in horizontal and 
vertical components, we arrive at the following model expression for the force f„ applied 
by a single filament pair onto the cell body center x^: 



l f A (hor) 



a 



(ver) 



2 

/(Ay 



,V f *+ f i) 

(cos fa + cos (pj)dij + a(sin fa + sin 4>j)djj 



(43) 



where a > is an additional adhesion or friction parameter. By relying on the pairing 



filament density p(9) in section 3.2 we obtain an integral expression for the total pair 



interaction force applied by cell i onto cell j: 



r (int) 
ij 



Rij / d# p{6) ■ 



(44) 



where the trigonometric relations between 



to be extracted from equations (30) 
determined by the relations 

^(hor) _ _ 



(33). 



^(hor) 



4>j and the parameterization angle 9 have 
Conversely, the force of cell j onto i is 



,i(ver) 



F, 



(ver) 



(45) 



The emerging cell pair interaction force as described by equation ( 44 ) is shown in figures 



12 and 13 A natural and maximal cut-off distance for the force is given by the finiteness 
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Figure 12: Modulus i*^- of pair interaction force F^ nt \ see equation (44). In the empty 



regions of the plot, is not defined. There, £>R i0 (xj) n Bj^ (xj) extends 

beyond as bounded by 6t, or the cells are not in contact at all. 
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Figure 13: Modulus of pair interaction force F^j depending on (a) 5ij and (b) n. The 
global view depending on both 5ij and the logarithmic ratio rj/r^ was pre- 



sented in figure 12 
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of the Voronoi tessellation, whereby neighboring is only possible for sufficiently small cell 
center distances 5ij < 5^ up ^ + r j + rj, i.e. Z3R i0 (xj) n B^ (xj) / 0. Once two previously 
isolated cells come close enough for contact, there is a strong tendency to attach, which 
facilitates multicellular tissue formation. The interaction force is attractive until the 
cell distance 5ij reaches <^ nt ^ = A cr j t • (>fj up \ where F^ nt ^ vanishes. Finally, if dij drops 

below F^ nt ^ becomes repulsive and therefore hinders tissue collapse at distances 

approaching S^J 1 = A mm • <5^ up . Note that with Aij > A m i n the lower bound from 



inequality (28) on the homogeneity of cell radii due to fixed V^max can be relaxed to 



1 < 7 ; ^ ~ V 2 < Qnb- (46) 

Correspondingly, P max can be increased for given cell homogeneity Q or Q n b- For exam- 



ple, the constraint (46) yields Q n b = 6.25 for y/V max = 3, or r m \ n > 0.73 • r max for each 
cell pair. In fact, the actual distances Ajj in a tissue will be larger than Amin 5 effectively 



relaxing (46) even further. 



3.4 Locomotion force at the free boundary 

In addition to the dynamics induced by pair interaction forces, cells at the tissue margin 
may migrate into open space. The locomotion force causing such a migration is due 
to lamellipodial protrusion and retraction, which is unhindered only at the free cell 
boundary Tjo- In a similar manner as before, we assume that this locomotion or free 
boundary force onto the cell body B ri (xj) is determined by connecting radial filament 
bundles as indicated in figure [9} The filament density of cell i along its free boundary 
r i0 is given by 

P ri P ( A1\ 

K i0 V ' max 

and thus independent of rj. In this way, the locomotive force of a cell i reads as 

F?° c) = /ioc J d Sl Pio&iO&i), (48) 

with arc length Si = Rio<pi and Rjo(^i) = (cos fa, sin fa). Moreover, in order to heuristi- 
cally account for ubiquitous perturbations due to lamellipodial fluctuations or possible 
signals, we implement stochastic force increments at the tissue margin 

dFf t} =b [ dB t>Si . (49) 



Here we assume a uniform and anisotropic vector noise B( jSi defining a spatio-temporal 
Brownian sheet in arc length and time coordinates with independent Gaussian incre- 
ments satisfying Var(dB( Si ) = dsi • dt. For each time t, stochastic integration results in 
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a simple weighted Gaussian noise term with random increments dWt 



dlf ) = b 0V ^\dW t = bo^^ptt, (50) 

where £t is a vector of Gaussian random numbers, which is chosen independently for 
every time step in a corresponding numerical realization of the stochastic process. 

3.5 Drag forces 

Apart from interaction and free boundary forces, the cell is subject to drag forces F^ drag ^ 
slowing down its movement. Such drag forces are generally functions of the cell body 
velocity Xj = Vj. Here we assume the simplest dependency of a linear force- velocity 
relation 

F f rag ) = - 7l v 4 , (51) 

with drag coefficient ji = 7(rj). Arising from friction with the substratum, could 
depend on the area of the cell body, e.g. 7(7*$) oc r|, however, for simplicity, here we take 
7i = 7 independent of cell body sizes. 

3.6 Dynamics of cell movement 

The previously described, active and anisotropic forces F^ nt \ f| 1oc ^ arising from the actin 
filament network act onto the cell center Xj causing translocation of the cell. However, 



friction, see equation (51 ), is considered to be dominating and inertia terms are neglected 
[2T1 l37| [25] . so that the emerging deterministic overdamped Newtonian equations of 
motion read as 

1 I —i(loc) V - ^ m(int) \ F,; 



F Uoc) + ^ F (mtM = .^ (52) 
v j neighbor ^ 

Moreover, any change of the translocation direction as well as adjustment of speed to 



the pseudo-steady state as given by the previous equation (52) requires some (mean) 
time Tj for restructuring and reinforcing the anisotropic actin network. The simplest 
way to model this adjustment process is by a linear stochastic filter of first order for 
the velocity [5|. Together with equation (49) this results in the stochastic differential 
equation (SDE) system 

dvj = ( — - Vj ] dt + bi^\Tio\dW t , dx t = v^dt, (53) 

-ti\7i / 

with bi = bo/^fi. Similarly as the friction 7^, also the mean adjustment time Tj could 
have some dependence on rj, however, here we restrict ourselves to the case of cells with 
homogeneous activity time scale V2 : Tj = T. 

For each time i, the forces ( 44|48|50|51 ) can be computed explicitly from the Voronoi 



tessellation of the generating cell bodies {i3 ri (xj) : i = 1 . . . N} using a spatial discretiza- 
tion of {Tij} in the parameterizing angle 9. Next, the velocities {v^} and positions 
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{xj} of the cell centers are updated according to both equations (53) in an explicit 
Euler-Maruyama step [28 . Finally, the Voronoi tessellation is computed anew from the 
updated generators {£> ri (xj)}. 

Higher order stochastic integration schemes were not applied, since such procedures 
necessitate the distribution of both forces and perturbations onto the powers of a Taylor 
expansion. In the general case, the involved derivatives of cell-cell contacts {T^} and cell 
margins {Tio} cannot be computed easily a priori. In particular, the change of a contact 
Tij may depend on the behavior of several distinct nearby cells k ^ Altogether, 
here we use the versatile basic method for integrating the equations of motion, because 
it is applicable regardless of the structure of the underlying SDE system. 



4 Cell pair contacts 

From the force plots in figure 
cell pair exhibits a sharp onse 
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it is clear that the interaction force ivj between a 
when two formerly dissociated cells come into contact. 



For any such cell pair in contact and each cell-cell body distance 5 = 5i 



equation (40), there is a unique pair of maximal contact angles (f>io,^)jo. 
according to equations (21) and ( |34| with r\ = Tj/ti and figure [7] the relation 



Therefore 



sec 



R j0 ■ sm((f>j (5)) = R i0 ■ srn(^ i0 (^))- 
holds and the free boundaries of both cells k = i,j are characterized by 



r fcl 



An 



R k o ■ (cos 4> k , sine 



> 



'/,-(> 



Thus, by solving the integral in equation (48) and regarding relation (47) we obtain for 



the deterministic part of the locomotive forces 



F (ioc) _ 2 f\ oc p n ■ s in( 
= 2/iocpr^ -sin( 



Mac) 



again using the decomposition in horizontal and vertical components. This means, that 
under our model conditions the mean locomotive forces of the two cells are exactly 
opposite, inde pen dent of their size. Since the interaction forces have the same property 



(see equation ( 45 ) with F^J = 0) , we conclude that the sum of the deterministic driving 



forces onto the cell pair vanishes, Fj + Fj 
onto cell j is computed as 



0. From equation (39) the interaction force 



■CI 



(5))d 



(54) 



with C((t>) 
In figure 
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Jq dtp (cos tp + \J\ — r^sin 2 tp) /^(tp) 1 / 2 , and k v as defined in equation (35). 



the resulting scalar horizontal force Fj 



^(int) +F 0oc) . g plQtted 

as a function 
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of cell body distance 5 and the locomotion force parameter /i oc , revealing the emergence 
of a stable deterministic contact equilibrium -F^det = for lower values of /i oc < 105 pN. 
On the other hand, for larger locomotion parameters, no such equilibrium exists and the 
cell distance always increases until the pair separates at 5 = S^ up ^ (= 10/im in figure 



14) 



In order to study the full stochastic contact and segregation dynamics described by 
the SDE system (53) with noise amplitudes bk = &o \/Gk{(t>ko{$))/jk, Gk{<t>) = \^ko\ = 



2\/7 5 maxT/c(7i" — <t>), we set all vertical noise components to zero, for simplicity. Then due 
to Fi + Fj = the dynamics is determined by differential equations for the cell overlap 
z and the difference u in horizontal cell velocities: 

z = S^ p) -5>0, u = vf° r) - v^ or) . (55) 

Proposition 2. For cell pair dynamics restricted to the horizontal connection line the 



stochastic ODE system in equation (53) transforms into a nonlinear Ornstein- Uhlenbeck 



system for the overlap z and its temporal change u, namely 



dz = udt (56) 
du = [f{z)-u^ + b{z)-dW t (57) 



Defining the mean drag coefficient '■= + and C{4>) as in equation (54), 

we have 

F(z) = ^(/(z)C(<A)-/ioc-sinA /(g) = / int log( W ~* Y (58) 

Ttj ^ ' ^ ^max Z c ' 

^)^( G f + G f), (59) 
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with z n 



r( ru P) 



5g up) • (1 - A r 



and z r 



? (rup) 



T/ie relation between 



5 and (j) = (j>jo(6) can be written as 



z = R j0 [l-Jl-sm 2 <t> \ +R i0 1- a/1-77 2 sin 



(60) 



Note that the overlap z is a monotone function of sin eft, which is proportional to the 
vertical extension of the overlap region as spanned by the contact arc Tjj (shaded in 
figure [7]) . 

4.1 Asymptotic stochastic differential equations 

Disruption of a connected pair occurs as z —> 0, so that an expansion at zero of all terms 



in Proposition 2 is justified. First, from equation (60) we derive the asymptotic relation 



z = R j0 ^-H s in 2 <t)( 1 + C(sin 2 



for z > so that 



sin i 



(1 + rj)R 



■jo 



Thus, the locomotion term of the force F(z) in equation (58 ) has a singularity at zero like 
\fz. The same holds for the interaction term; as a surprise, the corresponding integral 
can be expanded in (f> independent of the ratio r\ = Tj/ri'. 



C((p) = j dtp (2-^ 2 + e>(i/)) =2 -sin <£(l + 



0( 



The conclusion is that the horizontal pair force can be approximated as 
F{z) = 2^1 smJ2f(z) (l + 0{z A )) - / loc 

where only the prefactor depends on the cell body sizes. The deterministic equilibrium 
overlap z* > 0, as the zero of F, is approximately determined by the force equality 



/(**) 



loc 



Thus, by using the definition of / in equation (59), we define the contact parameter 

/loc 



A := z 



max V^niax 



z c ) exp 



2/ 



int 



(61) 



Let us consider the deterministic overdamped dynamics z = F{z) obtained from the 



stochastic ODE system (56, 57) in the limit T — > 0. With -fi = jj = 7 we prove the 
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Proposition 3 (Asymptotic contact and disrupture dynamics). Given a pair of cells 
with body radii r^ and rj, with small contact parameter |A|. Then there exists a unique 
stable equilibrium > if and only if X > 0, namely z* = A(l + 0(|A| 2 )) > 0. Moreover, 
in the limit T — ► 0, (3q = bo/^/T = const., and the corresponding stochastic differential 
equation (SDE) can be approximately written as 



dz = /V^fln r)dt + P(z)dW t . (62) 



Here z + = max(0, z) and 



G 
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with Gij = (8^V max ■ nrj/in + r,-)) . 



; ; / int , (3{z) = po \2^/l\ ^ ^{{r ^ + rj )ir - G l3 ■ Jz 



1/2 



Note that the log term in equation ( |62| ) can be approximated by (A — z + 0(z 2 )) /z max , 
and z* = A is a measure of the mean cell-cell overlap. 

As soon as the contact between cells is lost, z < 0, their cell center distance dij = 
VPm&x( r i + r j) — z would perform a pure Brownian motion for T = or, for T > 0, a 
persistent random walk with the standard recurrence probability to hit the touching state 
z = again. However, for situations of tissue cells moving on two-dimensional substrates, 
the production of adhesive fibers (as fibronectin) or remnants of plasma membrane plus 
adhesion molecules in the wake of a migrating cell would induce a positive bias of the 
locomotion force towards the other cell |27j . which could be assumed as proportional to 
the cell boundary distance —z, at least for small distances. Therefore, and for the aim 



of exploring the resulting stationary process, instead of equations (62) we consider the 
extended SDE model 

dz = F{z)dt + (3{z)dW t (63) 



with 

F(z) 



H\fz In I J for z > 

\ 2max Ay 

— VZ for Z < 



where we introduce an additional bias parameter v > describing an indirect attraction 
between separated cells. 
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4.2 Analysis of the stationary contact problem 

By solving the stationary Kolmogorov forward equation, we compute the approximate 
stationary probability distribution for the overlap z as 



z 

p(z) =p exp (J dz 
o 



a 3 /2 , A i 2 \ p (64) 
z2 -AH z2 z for z > 



V^max ^3 2x 5 

\z\ for z < 

2^ 

with a unique normalization factor po an d additional lumped parameters that arise from 
expanding the singular noise variance (3 2 {z) = 2-0(1 — sj z(t)/x + @(z 2 (t)): 

1p =\/^ , max/?o 7r ( r « + r j) 2 



,2 



1 1 



In figure 15 [b) and (e), the probability distribution p(z) according to equation (64) is 
plotted for two special parameter sets together with the force function F(z), whereas 
in figure |l5[c) and (f) histograms for the corresponding numerical realizations of the 



stochastic differential equation (63) are shown. For the contact parameter A = 0.25 > 0, 
see definition 61 , there appears a skew-shaped modified Gauss distribution around the 
unique center z* = A with a certain probability P Q fj = J ^ dzp(z) rj 0.15 for the cells 
to be separated. In contrast, for A = —0.15 < 0, the standard Gauss distribution with 
variance ip/v for positive separation distances — z is cut off by a decreasing distribution 
of the positive overlap z with mean value Z con t = dz zp{z) ~ 0.08 <C z* and a certain 
contact probability P C ont = 1 — -Poff ~ 0.32. 



Figure 16 shows the plot of the separation probability P Q g over the contact and noise 
parameters A and (3q together with the contour curve for the critical value A0.05 = A(/3o), 
so that the separation probability P Q g is less than 5% for contact parameters A > Ao.o5) 
i.e. for 

2/int ^ , ( A cr i t — A 
> log 



/loc ' Vl-A min -A .o5/4 mP) - 
For the two cases of contact parameters A > and A < 0, the corresponding time series 
of the overlap distance z as plotted in figure [l5|a) and (d), respectively, reveal a clearly 
distinct temporal separation and contact behavior. 

In the first case, with lower locomotion force parameter /i oc relative to the interaction 
parameter /i nt , the contact state persists for longer time intervals. Thereby, the overlap 
fluctuates around the equilibrium value = A, with mean contact duration T CO nt ~ 
exp(2.8) min « 16min. This proper contact time is evaluated from the maximum of the 
right-side r-hump in the logarithmic histogram of figure [l7|^a). The contact states are 
interrupted by periods of faster flickering, as can be seen from the change between on 
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Figure 15: Analytical distribution (bold/blue in b,e), the corresponding simulation his- 
tograms (c,f), and the time series (a,d) of the overlap width z. This quantity 



emerges from the dynamics of equation (63) for contact parameter A > 
(a-c) and A < (d-f), respectively, with force F{z) (grey/magenta in b,e). 
Time t is given in units of min, overlap z in units of //m, and force F in nN. 
Parameters used are (3q = 0.1, tp = (3$, z max = 0.1 and z c = 0.5. 
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Figure 16: Separation probability P g = f° oQ dzp(z) according to (64 1 plotted over the 



parameter plane of contact (A) and noise (/3q). In this plane the contour curve 
PoS = 0.05 is drawn in blue. The other parameters are the same as in figure 



15 More details are given in the text. 




log (contact time 



log (separation time) 



Figure 17: Logarithmically scaled duration histogram of contact intervals (a) and sep- 
aration intervals (b) for a longer time series (3000 min) of the situation in 
figure [l5T a). 



30 



and off with a mean duration of Tfli c k er ~ exp(— 1.5) ~ 15s, which is visible at the joint 
lower maximum in both figures [l7|a) and (b). Otherwise, an intermediate separation of 
the order of 1-5 min occurs. 

Yet for the second case, A < 0, with relatively higher locomotion force parameter /i oc , 



we observe in figure 15 'd) dominating periods of flickering around the steady state of 
touching = 0, which now is a stable deterministic equilibrium, even super-stable for 
z > with convergence z{t) ~ (i* — i) 2 / 3 in finite time. Longer periods of separation 
are also evident but their distribution does not differ much from the situation for A > 
(similar to the histogram in figure [TTJ^b), not shown here). 



5 Tissue simulations 

After these analytic considerations in the case of cell pair formation we return to the 



full equations of motion (53). Since both time and length scale of cell motility processes 
are well known, the only remaining free scaling figure is the magnitude of cell forces. In 
accordance to [3], here we assume that a typical bundle of several actin filaments can 
exert a force of approximately 10 pN. A single cell can, with the overall filament density 
parameter p and the force prefactors /ioc ; /int as in table [TJ reach an effective traction 
of 0(1000 pN) from a force as given by equation (|44[). The drag coefficient 7 then 



v'Pmax — 3 


p = 9.55/ pm 


7 = 2.5 • 10 4 pNs/^m 


dt = 2s 
T = 120 s 


A min = 0.1 
A crit = 0.2 . . . 0.7 


/int = 60pN 
/toe = 10... 20 pN 


6 -(7-Vax) 1/4 = 8.31plN/V/mis 


a = 0. . .0.17 



Table 1: Simulation and model parameters as described previously. Unless indicated 
otherwise, all simulations have been performed with this parameter set. 

naturally follows from experimentally observed cell velocities |50l [20) [33] . As explained 
before, T is the persistence time of the intracellular cytoskeletal reorganization, and ? max 
determines the relative size of the lamella region around the cell body ri (xj). Moreover, 
the scaled interaction distances A m i n , A cr i t as defined in equation (42 ) determine the sign 



and scaling of the cell pair interaction force F^*; . The relative strength of the vertical 



component of F?- in equation (43) is given by a. Finally, the stochastic perturbation 



parameter bo in equation (50) contains a factor (7-niax)~ m order to obtain the same 
amount of perturbation for cells with equal body radii rj. Since we look for robust 
features in the simulations, &o was chosen fairly high. 



5.1 Emergence of tissue shape and multiple stable states 



Consider a simple proto-tissue of seven cells as shown in figure 18 'start'. Using the 
parameters /i oc = 20pN,a = 0,A cr i t = 0.25, a series of 1000 simulations has been 
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(b) 



(c) 



(d) 




(b) (c) (d) 
conformation 



Figure 18: Different tissue conformations (a-d) evolving from the configuration 'start' 
after a simulation time of 8h; here r max = 2.0 fj,m (cell 7), r m i n = 1.0 /im (cell 
2), and parameters /i oc = 20 pN, a = 0, A cr it = 0.25. The percentage of 
occurrence of a particular conformation then was computed, and the error 
bars were obtained by a simple bootstrap method. Other features are further 
explained in the text. 
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performed. After a simulation time of 8h, the emerging tissue conformations as distin- 
guished by Delaunay network topology have been recorded. In the course of these 8h, 
significant changes appear within the tissue, and apparently several equilibrium confor- 
mations emerge. For the four most prevalent conformations the percentage of occurrence 
is displayed in figure 18 One observes two rather globular shapes (a), (d), where either 
the big cell 1 or the two small cells 2, 3 are engulfed by the others, respectively. Further- 
more, there are two more elongated shapes (b), (c), where only the single small cell 2 is 
completely surrounded by other cells. Being distinguished by topology, (b) and (c) are 
in fact quite close in shape, despite of their rotational variation. 

From the high occurrence of the topological conformations (a), (b) one might conclude 
that these two conformations are the most stable ones. Thus, and to clarify the interrela- 
tions between the conformations (a-d), we investigate (a) and (b) in longer simulations. 
To this end, by starting from the configurations (a) and (b) (see figure 18), both tissues 
have been evolved for 40 additional hours of simulation time. In order to characterize 
the shape of tissue with respect to global and elongated shape, here we observe tissue 
size, i.e. the maximum diameter, and tissue circularity 

_ 2 a/ 7t A tiss < 

J2i l r »ol ~~ 

where Ati ss is the total area of the tissue. Note that in connected tissues = 1 would be 
attained for a purely circular globe. Other observables are possible but less indicative 
in this context. In figure 



19| the two time series (blue, red) soon reach the different 
conformations of figure 18 a) and (b) respectively, at times around t ~ 8h. While the 
circularity Q is only slightly different in the two cases, the tissue size is clearly higher 
for the elongated conformation (b). Apart from stochastic fluctuations and an initial 
equilibration phase for t < 1 h, both observables attain a constant value for time series 
(a). In contrast, for time series (b) there appear distinct states between t ~ 2.5 h 
and t ~ 24 h. Indeed these observations are reflected by the actual evolution of the 
tissue. While the topology of the tissue does not change after t = 2 h for time series 
(a) (see supplementary |mova . avi| ) , tissue (b) ( |movb . avi| ) goes through several different 
conformational states. At t = 13.5 h it attains the same topology as conformation (c), 



identifying (c) as a transient state (movc.avi). Afterwards, approximately at t ~ 18.5 h 



cell 7 establishes contact with cells 1,3, so that the tissue shape is similar to (c). Finally, 
shortly before t = 24 h, the tissue reaches its final conformation similar to (d) except for 
the order of the marginal cells. Conformation (d) emerges in a similar manner as (a), 
however instead of cells 1,5 initially cells 3,4 form a neighbor pair, quickly leading to 
the stable final arrangement in less than 0.5 h (see supplementary movd . avij ) . Moreover, 
the time series of conformation (b) in figure 19 suggests, that during t = 0.5 . . . 3h the 
tissue already attains a shape of similar compactness and stability as in figure (18^ a). 
Nevertheless, the Delaunay mesh (green lines) is not convex there (movb.avi), which 
explains this surprising instability. 

It appears, that the stability of a tissue is related to its globular shape. This is 
not a surprise, since the stochastic forces in equation (48) are defined only on free cell 
boundaries Tjo, and therefore act only on marginal cells. Thus, by minimizing the extent 
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Figure 19: Two time series (blue, red) for the tissue in figure 18 'start'. After 8h the 



two distinct topological conformations from figure [l8[a) (blue) and 18 [b) 
(red) have emerged. These topological conformations are characterized by 
distinct tissue size (right axis) and circularity 0, (left axis). While (a) is 
apparently a stable topological conformation that does not change even for 
strong stochastic perturbations, (b) relaxes into a topological conformation of 
globular shape via several intermediate steps, see supplementary 
movb . avil 



mova. avi 



of all Tio, a maximal circularity minimizes stochastic perturbations, which enhances the 
stability of the tissue. Additionally, for a tissue to change its topology, its cells have to 
overcome barriers as imposed by the other cells. For example, cell 3 has to displace cell 
1 in |movb.avi| at t ~ 3h in order to make contact with 2. Depending on the particular 
configuration, the severity of these barriers might range from prohibitive to practically 
non-existent. Influenced by the strength of the stochastic interactions, these barriers 
then determine the time scale of further relaxation to equilibrium. In this sense, the 
notion of equilibrium is directly related to an inherent time scale. According to the 
previous evolution of the tissue, there may be several stable topological conformations 
for a given time scale. 



5.2 Stability of tissue formation 

In order to explore the ramifications of piecewise spherical cells within our model frame- 
work, we study the influence of "P m ax and A cr i t on tissue formation. To this end, a 
simulation has been performed starting from an exemplary configuration as in figure 
[6] with A cr it = 0.3, \/V max = 3, a = 0.17 and f\ oc = 10 pN. After 8h, either A cr i t or 



VP m&x was modified to a nearby parameter position as indicated in figure 20 and the 
simulation was continued for further 8h. This procedure was repeated until the whole 
panel in |20| was filled with the final tissue configurations. 

For fixed {rj}, the overall size of the tissue is determined by both v^Pmaxi defining 
the free cell radius in units of r%, and A cr ; t , presetting the equilibrium cell-cell body 
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V^W 1.5 2.0 3.0 5.0 




Figure 20: Stability of tissue for various values of the parameters y/V max . and A cr i t , where 
/loc = 10 pN, and a = 0.17. When increased, both parameters lead to an 
increased tissue size. For sufficiently large A cr it, the tissue eventually disso- 
ciates. The extremal cell body radii are r m ; n = 0.9 //m and r max = 1.7 /xm in 
all simulations presented in this figure. 
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Figure 21: Tissue aggregation and dissociation depending on the relative lamella width 
W and the dimensionless cell overlap 2. The line Z ■ W = To ~ 0.24 ± 3% 
separates associating from dissociating tissues and was determined from a 
linear regression against Z = Tq/W. See also the explanations in the text. 



distance in units of 5^ up \ see equations (21), (42) and (40). Correspondingly, in figure 
20 the overall tissue extension increases from left to right and from top to bottom. Fur- 



thermore, for given A cr j t , tissues with higher vPmax exhibit a rather compact, almost 
quadratic shape. We speculate that this is due to spontaneous formation of distinct pro- 
trusions arising from stochastic perturbations and leading to an increase of locomotion at 
the corners. In contrast, tissues with lower V^max feature more irregular margins. Sim- 
ilarly, for given V'Pmax, larger values of A cr it yield more irregularity, most pronounced 
directly before dissociation of the tissue. 

Apparently, the emerging interaction forces are sufficiently strong for tissue aggre- 
gation only if there is enough space for the adaptation of neighboring cell lamellae. 
This space serves as a cushion for accommodating near-range repulsion from multiple 
neighbor cells and at the same time poses a resistance to stochastic perturbations by 



mid-range neighbor cell attraction, see figure [14} Otherwise the tissue dissociates, lead- 
ing to isolated cells exclusively driven by stochastic perturbations. In order to quantify 
these findings, consider the relative lamella width W and the dimensionless cell overlap 
Z defined by 

V^max 1 



< 1, 



Z = l- A crit < 1. 
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that are not 



Inserting the marginal values V^max and A cr ;t of those tissues in figure 
yet dissociated, one recognizes that the product Z ■ W =: Tq is approximately constant, 
with Tq ~ 0.24, see figure 21 This Tq can be identified as a threshold value guaranteeing 



tissue coherence under the condition 



Z ■ W > Tq, 



(65) 
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where Tq eventually depends on the other parameters, which were fixed here. This 
confirms that for tissue aggregation to occur, V^Vax and with it the relative lamella 
width W must be sufficiently large. 

On the other hand, we have established this result under the tissue homogeneity 



condition (46) guaranteeing starlikeness of cells, which now can be rewritten as 



^max-W = l- ? ^, (66) 

with i? max := 1 — A mm defining the maximal dimensionless overlap. Since W < 1 by 



construction, inequality (66) always holds for very high cell size homogeneities Q nD > 
1/A m j n . However, for lower Q n b there is an upper bound on W restricting the available 
space for the cell lamellae. If, in addition, starlikeness of cells is enforced for all possible 



neighborhood constellations, then Q n b has to be replaced by Q, see equation (24). In 



this way the relations ( 65 ) and ( 66 ) lead to the sufficient condition for tissue coherence 



1 A C rit "v/T'max (1 A mm )(l + r max /r mm ) 

From these estimates we finally conclude that for given model force parameters (A m i n , 
A cr it, /int, /loo bo) the formation of integer tissue aggregates with overall starlike cells 
is guaranteed within a certain finite range of the free cell size parameter v^maxj where 
the upper bound decreases with an increasing ratio r max /r m ; n of extremal cell body radii. 



Within the limits of inequality (67) the lamellae regions are wide enough to perform 
the necessary deformations by adapting to the surrounding neighbors through shape 
changes. Thus, nature's freedom in developing aggregating tissues may be constrained 
by a tradeoff between the relative size of cells with respect to their bodies (V^max) and 
the cell size heterogeneity (1/Q). 

6 Results and discussion 

In this article we have investigated the emergence of tissue aggregation using finite, 
generalized Voronoi neighborhoods as a basis for the description of cells within epithelial 
tissues. It was shown that the two-dimensional geometric structure observed in tissues, 
in particular the circular shape of contact arcs and the size distribution of cells can be 
captured by the experimentally accessible characteristics of cell body radii {rj} and the 



relative extension P max of the lamellae. The quotient method defined by equation (17) 
with unique weight factors Wi = ri implements the expected partition of influence regions 
and additionally suggests the definition of directed force densities on the contact and free 
boundaries. In contrast, the difference method is not appropriate, because the influence 
region of a cell growing in size becomes reduced within a cell pair. Nevertheless this 
method has been used previously for modelling three-dimensional tissue dynamics |35[ 
[37] . By using the particular additive weights Wi = ri, the specific geometrical properties 
( "orthocircles" ) allow for a dense regular triangulation of space between neighboring cell 
bodies even in the case of overlap, see [TO] for details. 
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In the asymptotic case of two cells in brief contact, i.e. 5y ~ S^ np \ the approximate 
equations of overdamped motion without persistence were treated analytically. Thereby, 
the contact singularity at the rupture points could be resolved and a stationary proba- 
bility distribution was found under the assumption that separated cells indirectly attract 
each other due to biased locomotion. Corresponding simulations confirm that in the case 
of dominating attraction, i.e. contact parameter A > 0, there is indeed a unique, stable 
equilibrium position facilitating cell-cell attachment after the first contact. 

The simulations of multicellular dynamics, cf. figure 20, reveal a surprising richness 
of tissue structure, with emerging shapes from elongated to globular and from compact 
quadratic to irregular and dispersed. Ranging from widely spread to tightly contracted, 
different lamella protrusions appear for varying "Pmax- Similarly, within a more or less 
compressed tissue as in figure [6j single cells attain varying forms from almost circular 
at the margins to polygonal in the interior. Moreover, distinct stresses throughout the 
tissue lead to both narrowly compressed or widely stretched cell shapes. 

Furthermore, due to the stochastic nature of neighbor constellation and cell-cell inter- 
actions, different stable tissue configurations may evolve from the same starting configu- 
ration. The underlying dynamical time scales vary over several orders of magnitude. This 
can be interpreted in terms of an incomplete balance of competing interactions within 
a system. In the course of relaxation to equilibrium, this is known to cause frustration, 
meaning that the system ends up in a local minimum of the involved (generalized) free 
energy [38| [M] [26]. Thus, the structure of the proposed anisotropic and active force 
interactions exhibits interesting non-trivial features, yet it is sufficiently elementary to 
allow for a rigorous treatment under certain further assumptions. 

Finally, the limits of aggregation due to a prohibitively small tissue coherence threshold 
Tq have been explored. This leads to our main result following from purely geometric 
arguments: In proposition [T] we have derived that relatively large lamellae extensions 
"Pmax are restricted by the cell size homogeneity Q n b when requiring starlikeness of cells. 
On biological grounds one might argue that not all observed cells are starlike (in the 
mathematical sense). In such a case, however, the exertion of forces onto neighboring 
cells is certainly hindered - especially in cell regions that cannot be reached by radial 
filaments. There, the necessary centro-radial support of the apical actin cortex may be 
weakened due to the necessary bending of filament bundles. Reversing the argument, 
for the cell to provide macroscopic stability within the tissue, it is beneficial to attain a 
star like shape. 

Several model refinements are of interest. On the cell biological side, cell division is 
a commonly observed phenomenon, having special consequences for pattern formation 
and self-organization during embryonic development. In this direction also the growth 
of cells, or even the complete cell cycle including necrosis or apoptosis could be con- 
sidered. On the mathematical side more elaborate stochastics, such as differentially 
modeled stochastic forces on the free boundaries as well as on the contact borders could 
be implemented. In particular, in the light of recent experimental results |40|, 129] . vary- 
ing filament orientations or binding and unbinding processes of linker molecules might 
be considered. Moreover, we emphasize that a three-dimensional generalization of the 
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multiplicatively weighted Voronoi tessellation is straight-forward, opening a wide poten- 
tial for applications, for example epidermal tissue organization in the intestinal crypt 
without artificially imposed constraints. Finally, explicit variables for cell polarization 
and reorganization of the cytoskeleton could possibly lead to refined dynamics closer 
resembling the behavior observed in vivo. 
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